Coherence and Spatial Resolution of Transport in Quantum Cascade Lasers 



o 
o 

c 

oo 



c 

o 
o 



> 
o 

O 

o 



X 



Andreas WackeiQ 
Mathematical Physics, Lund University, Box 118, 22100 Lund, Sweden 
(Dated: February 1, 2008, submitted to Proceedings of HCIS-15 (July 2007), physica status solidi (c)) 

The method of nonequilibrium Greens functions allows for a spatial and energetical resolution 
of the electron current in Quantum Cascade Lasers. While scattering does not change the spatial 
position of carriers, the entire spatial evolution of charge can be attributed to coherent transport 
by complex wave functions. We discuss the hierarchy of transport models and derive the density 
matrix equations as well as the hopping model starting from the nonequilibrium Greens functions 
approach. 

PACS numbers: 05.60.Gg,73.63.-b,73.21.Cd 



INTRODUCTION 

Since the first realization in 1994 [l[ Quantum Cascade 
Lasers (QCLs) have become an important tool for IR- 
spectroscopy. The extension towards the THz region 
has opened up possibilities for a variety of applications in 
rapid-but-precise hazardous chemical sensing, concealed 
weapon detection, non-invasive medical and biological di- 
agnostics, and high-speed telecommunications 0. 

The operation of QCLs is based on electronic transi- 
tions between different subbands within the conduction 
band of a semiconductor heterostructurc. Using a sophis- 
ticated sequence of wells and barriers, the electrons are 
guided into the upper laser level at the operating bias, 
thus creating population inversion for a pair of levels in 
the active region. By modifying the layer thicknesses, the 
transition energy can be varied in a large range and oper- 
ating lasers with wavelengths between 2.95/Ltm [ij and 217 
jUm (1.39 THz) [B| have been realized within the last year. 
This covers a large part of the electromagnetic spectrum 
from the near-infrared to the proximity of fast electrical 
circuits (albeit there is a gap at the Reststrahlenband) . 

Conventionally, QCLs are modeled by rate equation 
schemes either for the average electron densities in the 
subbands [fl M, |8fl, or the occupation of the individual 
states 0, fiol. TllL 12 . Il3j |. The latter ones arc often 
simulated with the Monte-Carlo technique, which sug- 
gests this denomination, albeit the term hopping trans- 
port [3, EH seems to be more appropriate for this type 
of models. Such simulations allowed for a continuous im- 
provement of device performance by optimizing the layer 
structure for an appropriate ratio of scattering matrix el- 
ements and resonance conditions. Hopping or rate equa- 
tion models can however not describe coherent effects, 
which are of some relevance for the tunneling transition 
between the injector into the upper laser level 

nana. 

Furthermore, the broadening of the gain transition can 
only be qualitatively estimated within such models. To 
overcome these limitations, a quantum transport model 
based on nonequilibrium Green functions (NEGF) was 
developed [H| • It was demonstrated that the microscopic 
current flow is due to coherent evolution of wave packets 



rather than the spatial translation by scattering transi- 
tions Here this idea is further elaborated with a par- 
ticular focus on the relation between the NEGF model, 
density matrix equations [2fj| . and the above mentioned 
hopping transport models. 

The paper is organized as follows: In Section 2, the 
different concepts for calculating a current in QCLs (or 
similar semiconductor heterostructures elements) are dis- 
cussed. A key result is that the entire current is carried 
by nondiagonal elements of the density matrix (coher- 
ences), as also discussed in Ref. [2(|. In Section 3, we 
present numerical examples for the different representa- 
tions of current. In the more technical sections 4 and 5 
it is shown how the density matrix equations, and the 
standard hopping models are derived by successive sim- 
plifications of the Greens function technique. 



MODELING THE CURRENT 

In planar semiconductor heterostructures, such as 
QCLs or supcrlattices, it is appropriate to use a set of 
normalized basis states -^=(^ a (z)e lk r which separate the 
behavior in growth direction (z) with quantum number 
a from the plane wave behavior (k) in the (x, y)-plane of 
total area A. The Hamilton operator is written as 



H= H° 

=P^t(T)P+ V <:( Z ) + e <f>( Z ) 



-H s 



(1) 



where V c (z) = V c (z + d) is the conduction band edge for 
our structure with period d. 4>{z) is the (self-consistent) 
electric potential satisfying 4>(z + d) = <f){z) — Fd, where 
F is the average electric field along the structure. Here 
it is important to note, that H° is diagonal in k due to 
the translational symmetry of the perfect QCL structure 
in the (a;, y)-plane. In contrast, impurities, phonons and 
possibly the presence of other electrons constitute scat- 
tering terms of the form a^ k aak' with k ^ k' in H sclLtt . 
Here a ak and a a k are the standard creation and anni- 
hilation operators in occupation number representation, 
respectively. In the following we use the Wannier basis 
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for our calculations, see [21|, which provides a periodic 
array of states satisfying ip a (z + d) = ip a i{z). The same 
property holds for Wannier-Stark (WS) states as well, 
which in addition diagonalize H° with energies E a (k). 

The current can be evaluated in two ways: The current 
density averaged over the entire sample reads 



2 (for spin)ie ^ 

a,/3,k 



(2) 



hv 



where p a p(k) = (a^ k a Q k) is the density matrix (here de- 
fined to be diagonal in k), Wp a = J2~< H p~< z ~fa - z/3yH° a , 
and V = NdA is the normalization volume of the sys- 
tem with N periods. Note that the second contribution 
of the Hamiltonian fTj) provides {[H scatt ,z]) = 0, as the 
operator H 8catt is only a function of f, but not of p, for 
all scattering processes typically considered flit ]. 
The local current density is given by 

J(=) =«/ dlds _L 



h d 



h d 



mA 



£(4k«/3k) 



(3) 



a/3k 



x \ ^(^)t^(z) + 



V»(*) 



where the expansion ^(f) = 5Zg k </'/3(-2 ; )e lk ' r a^/ £ /v / 3 for 
the field operators was used. Averaging over z and using 
the commutator relation 
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m(z) dz dz m(z) 



(4) 



provides directly Eq. @. 

It is important to notice, that the matrix W a /3 is anti- 
hermitian. Thus the diagonal elements vanish for a set of 
real basis functions and consequently the entire current 
is due to the nondiagonal elements of the density matrix 
p a p(k). The same holds for Eq. Q: If the wave functions 
ip a (z) are real, the diagonal elements of the density ma- 
trix do not provide any contribution to the current. Now 
both the Wannier and Wannier-Stark basis functions can 
be chosen real and therefore in both cases the current is 
entirely being carried by the non-diagonal elements of the 
density matrix p Q) g(k). 



by 



Working with NEGF [22j, the density matrix is given 



(5) 



Thus, the correlation functions G < (E) can be viewed as 
the energy-resolved density matrix and Eq. ([3]) can be 
generalized to the energy-resolved current density 



a/3k 



H > f \ 



(G) 



This equation becomes of particular interest, if one 
considers a special basis set of states ^^(E, z), which 
diagonalize Gp a (k, E)/(2m) with the real (and positive) 
eigenvalues f n k{E). Then the current as well as the den- 
sity is represented by an incoherent superposition of com- 
plex wave functions ^f n k(E, z) at each energy. The fact 
that these wave functions carry the entire current man- 
ifests the coherent nature of current evolution in QCLs 
as well as related structures such as superlatticcs. 



NUMERICAL EXAMPLES 

Let's consider the THz-QCL from Kumar et a/.[23|. 
which operates above 77 K. The current-voltage charac- 
teristic evaluated via Eq. (|2|) is shown in Fig. [TJa) and 
good quantitative agreement with the experimental data 
is found. (Details of the calculation are given in [3].) 
In Fig. [Ijb) this current (dashed line) is compared with 
the local current density (full line) evaluated via Eq. ([3]). 
While current continuity requires a constant J{z) in the 
stationary case, the evaluated local current exhibits spa- 
tial oscillations. The amplitude of these oscillations de- 
creases with the number of Wannier states per period 
employed in the calculations. This suggests that this ar- 
tificial effect is due to the lack of completeness if only 
a finite number of basis states is taken into account for. 
The result from Eq. ([2|) corresponds to the spatial average 
and is far less sensitive to the number of states employed. 
This shows that the current evaluated by Eq. (PJ]) has to 
be taken with care and Eq. ([2]) is preferable. 

In the upper panel of Fig. [2] the energetically resolved 
current density from Eq. ((6J is displayed. At each en- 
ergy one observes a current flow in z direction due to the 
presence of nondiagonal elements in G^ a (k, E) . In order 
to satisfy the continuity of current, scattering transitions 
transfer particles between different energies, where the 
coherent evolution of the current continues. In addition, 
there are also elastic scattering events, where the cur- 
rent continues at the same energy, but with a different 
parallel momentum k, which are visible in corresponding 
k-resolved plots. 

For comparison the energetically resolved particle den- 
sity (see also [Hj]) 

n(E,z) = j^±G$ a (k t E)<p* a (z)<p fi (z) (7) 

a/3k " 
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FIG. 1: (a) Current-voltage characteristic for the THz-QCL 
of ref. [23J . In the calculations (full line) the bias was taken as 
177 times the voltage drop per period and the area as A — 54- 
10~ 5 cm 2 . The cross marks the operation point at a bias drop 
of 62 mV per period. The experimental data (dashed line) 
are by courtesy of S. Kumar, (b) Spatially resolved current 
density from Eq. © for a calculation with 5 and 8 Wannier 
states per period, respectively (full lines). The dashed lines 
give the result from Eq. ((2} for comparison. 



is shown in the lower panel of Fig. [2] It is intriguing 
to sec, that neither the density nor the current profile 
follow the spatial profile of the WS-states. Furthermore 
note, that the states 1 and 2 resemble the binding and 
anti-binding combination of two more localized states. 

In Fig. [3] the eigenvalues of G^ a (k, E)/(2iri) are dis- 
played as a function of energy. One can identify distinct 
peaks indicating the presence of broadened quasiparticle 
states. These can be attributed to specific branches in 
the eigenvalue spectrum, which however mix with each 
other at crossing points. The essential structure for k = 
is repeated for finite k-values with a shift in energy by 
E k = ti 2 k 2 /2m. For E k = 10 meV the width of the 
peaks is larger as more scattering states are present than 
for Ek = meV. 

In Fig. S{a,b) the wave functions corresponding to the 
three largest eigenvalue peaks are shown. They describe 
the spatial structure of both the electron density and cur- 
rent displayed in Fig. [51 Thus they give a better descrip- 
tion of the ongoing behavior than Wannier or Wannicr- 
Stark states. In many cases these states are essentially 
unchanged if one follows a single branch of eigenvalues in 
Fig. [3J E.g., the states 1 and 4 differ only very slightly, 
see Fig. [He, d). However a strong energy dependence can 
occur due to mixing effects if different branches of eigen- 
values come close to each other or even cross. This can be 
seen in the sequence for states 5, 6, and 7. The states for 
finite k are related to the corresponding states at k = 
at the lower energy E — Ek- E.g., state 8 corresponds 
to states 5,6 (which are about 10 meV lower in energy 
E), albeit the mixing between different branches makes 
a detailed comparison difficult. 
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FIG. 2: Upper panel panel: Spatially and energetically re- 
solved current density evaluated by Eq. ([6]). The WS states 
<Pa(z) corresponding to the upper (1') and lower (5) laser level 
are depicted for orientation. The vertical array marks a rep- 
resentative scattering transitions. Lower panel: Spatially and 
energetically resolved particle density evaluated by Eq. ([7|. 
The lowest five WS states tp a (z) are displayed. 



DENSITY-MATRIX EQUATIONS 

Now we want to study the relation between the dif- 
ferent approaches. In the NEGF approach, the results 
are determined by the lesser Greens function. For the 
stationary state, where 

G<(t,t')= f ^G < (E)e-^ t - t "> E ^ h 7 

Eq. (5.4) of Ref. [H] provides us with 

^G< 7 (£,k)i/yk) - H° a7 (k)G<p(E,k) 

7 

= £ [^(E,k)G< (E,k) - G™*(E,k)X< (E,k) (8) 

7 

+ E< (E, k)G^(E, k) - G<(E, k)Z»f(E, k) 
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FIG. 3: Eigenvalues of the matrix G5 a (k, E)/(2iri) as a func- 
tion of energy for two different values of k. The states corre- 
sponding to the eigenvalues denoted by circles are displayed 
in Fig. H 
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FIG. 4: Wave functions >]/nk(i5, z ) which diagonalize 
Gp a (\i, E) / (2ni) for different energies and k corresponding 
to the eigenvalues denoted by the circles in Fig. [3] 



The self-energies are evaluated in the self-consistent 
Born-approximation providing 

S< /3 (S,k) = ^F a7 (q)^(-q) 
7<5q 

x[nqG< s (E - fiw^k- q) 

+ (n_ q + 1)G< (E + hw^ k - q)] 

(£,k) = ^V Q7 (q)%(-q) 



-,rct/adv 



(9) 



x [(n q + 



^i)c; c ; /adv (£;-^ q ,k-q) 

n _ qG ret/adv (£; . 



fio;_q,k-q)] 



For illustrative purpose only phonon scattering with a 
single lateral mode q is taken into account here and the 
nondcgcncrate case is considered (otherwise additional 
terms with G < appear in j] adv A ct ). However, neither 
of these simplifications was performed in the numerical 
examples discussed above. 

Neglecting any broadening effects, the full Greens func- 
tions can be approximated by the bare Greens functions 



^<rct / adv 



{E, k) stf^ 



1 



E - E p (k)±i0+ 



(10) 



G<JE, k) K2mp af} 8{E - E a/3 (k)) 



A key issue is that we allow for a nondiagonal density 
matrix, which makes it difficult to address a specific en- 
ergy E a p(k) to the respective 5-function. A first guess is 
that E a p(k) is somehow related to E a (k) and/or Ep(k). 

Now Eq. ((9]) is inserted into Eq. ([8]) and subsequently, 
the approximations (|10|) arc inserted in the right-hand 
side. Integrating over E and dividing by 27ri, provides 



5> Q7 (k)# 7 Vk) - ff° 7 (k)p 7/3 (k) 

7 



7<5q 



-E 70 (k) - E s (k - q) + fi^-q + 

nqVgy(<l)Py6{k - q)Vtf^(-q) 
£ 7<5 (k - q) - E a (k) + huj cl + 10+ 

ngV a s{q)psy(k - q)Vyp{-q) 
E Sl (k - q) - E, 3 (k) + Hujq - i0+ 
(k)y 75 (q)%(-q) 



iOH 



(11) 



J E a7 (k) -E s (k-q) + huj-ct - i0+ 
terms with n a — > n_„ + 1 and fku n — > —Huj. 



Setting £ 7/3 (k) = £ 7 (k), ^(k - q) = £ 7 (k - q), 
£7<5 7 (k — q) = E 1 (k — q), and _E Q7 (k) = -E 7 (k) in the 
subsequent lines on the right-hand side, one obtains pre- 
cisely the density-matrix kinetics of Sec IID of [2(| in the 
so called complete collision limit. In this kinetics, the 



left-hand side has the additional term ih 



dp Q)3 (k) 
dt ' 



which 



however vanishes in the stationary case considered here. 
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In the density matrix equations, the choice of E a p(k), 
which accompanies the density matrix p a p(k) on the 
right-hand side, can be related to the way, the Markov 
limit is performed. Here different choices have been 
suggested prj 27 1, which is an issue of ongoing debate. 



However, as shown below, the nondiagonal density ma- 
trices are small unless |^g(k) — £? Q (k)| < T. If the prop- 
erties of the system are constant on this energy scale, e.g., 
the temperature is larger than T/ks, the specific choice of 
E a p(k) within the energy interval [^(k), i? Q (k)] is not of 
central relevance. Thus, the results for different choices 
should not differ dramatically as observed in [23] ■ In the 
opposite case of small temperature (< T/ks), broadening 
effects become of importance, which renders the density 
matrix approach questionable anyway. 



HOPPING MODEL 

The ambiguity of choosing E a p(k) vanishes, if we as- 
sume that the diagonal density matrices ppp(k) = fp(k) 
dominate the scattering terms which constitute the right- 
hand side of Eq. (jTTJ) . This makes particular sense, if the 
states are chosen as the eigenstates of H°. For a = (3 we 
find 

= -2m^7i_ q |^(q)| 2 

<5q 

x S(E a (k) - E 5 (k- q) + fiw_ q )/ a (k) 
+ 2^i]Tn q |K> 7 (q)| 2 ( 12 ) 

x S(E 7 (k - q) - E a (k) + ^ q )/ 7 (k - q) 
+terms with n q — > n_ q + 1 and fiw q — > — ?kj_ q 

This is just the difference of out-scattering and in- 
scattering transition rates for the state (a,k), where the 
scattering rates are evaluated by Fermi's golden rule. 
This defines the hopping model 11411 which has been fre- 
quently applied to QCLs @, [l(], EL H, 13j - It is usually 
solved by the Monte Carlo technique and provides the 
stationary occupations / Q (k). 

For a 7^ (3 the left-hand side of Eq. (TiTj) provides the 
term [Ep(k) — E a (k)]p a p in the eigenstate basis. Again 
the right-hand side has the magnitude of T x 0{f a }, 
where T/H is the magnitude of the scattering rate for 
a single level. Therefore does the assumption, that the 
diagonal elements dominate the density matrix, become 
questionable if a pair of levels satisfies | Ep (k) — E a (k) | < 
r which is typical for level crossings, see also the discus- 
sion in [l7j |. 

The evaluation of the current is a subtle issue, as the 
current is entirely contained in the nondiagonal density 
matrices as discussed above. Now Eq. gives in the 



eigenstate basis: 

J = ~W 2 ^{[Ep(k)-E a (k)}p a p} (13) 

where the antisymmetry of W a p was used. Now [Ep (k) — 
E a (k))p a p is precisely the left-hand side of Eq. (jlip and 
restricting to the dominating diagonal density matrices 
on the right-hand side we obtain 



J = 



27re x - 



Zpc 



\^2n- n V aS (q)V S p(-q) 

x S(E p (k)-E s (k-q)- 
- ^n q V r Q7 (q)V 7/3 (-q) 



^- q )//3(k) 



(14) 



X S(£ 7 (k - q) - E a (k) + fiw q )/ 7 (k - q) 



+terms with n q — > n_ q + 1 and hujq — > —hui-q 

where we used that the lower two lines are the complex 
anti-conjugate of the upper two lines in the right-hand 
side of Eq. (fTTj) after exchanging the indices a and (3. 
Now the completeness of the states <p a (z) in the z-part 
of the Hilbert space provides the relation 



zp a V a s(q) 



Ei 

/3 



\zV(z,q)\S)-ZppVps(q)} 



to be used in the first summand of Eq. (|T4|) . In addition 
the running index d is replaced by a. Correspondingly, 



E Vyp(-q)zp a 



is used in the second summand with the replacements 
7 — ► (3 as well as k q — * k and q — ► q. These 
operations result in 



J= W E n -v\ v Pct(.q)\ 2 {z aa - zpp) 

a/3,k,q 

x 5(E p (k) - £ 7 (k - q) + ^_ q )/ /3 (k) 



(15) 



-terms with 



- q + 1 and fiWq — ► —huj-q 



which is the standard expression for hopping currents. It 
can be interpreted as the sum of scattering transitions 
from (3 to a, which change the mean location of the elec- 
tron from zpp to z aa . This is however not the underlying 
physics, as scattering does not directly change the par- 
ticle position. In contrast the entire current is carried 
by the polarizations p a p and Eq. (|15[) is nothing but an 
approximation for these coherences. 
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CONCLUSION 

The transport in QCLs and similar structure such as 
superlattice is entirely due to coherences, i.e. nondiago- 
nal elements in the density matrix p a p (k) if a set of real 
basis functions is chosen. The use of NEGF allows for a 
spatially and energetically resolved visualization of these 
coherent transport properties. Neglecting the energetic 
broadening of the states, the density matrix equations 
can be derived from NEGF theory. However, the ener- 
getic location of the nondiagonal elements is only poorly 
defined in this reduction scheme. If the level differences 
are larger than the scattering induced broadening F, the 
nondiagonal elements of the density matrix are small and 
can be approximated by differences in level occupation. 
In this way the frequently used hopping model for the 
current appears. This model suggests the interpreta- 
tion that the spatial position of the particles is directly 
changed by the individual scattering processes. However 
one has to keep in mind, that conventional scattering 
processes do not change the position of carrier, but only 
induce coherences which subsequently drive the current. 

The author thanks F. Banit, A. Knorr, S.-C. Lee, R. 
Nelander, M.F. Pereira, C. Weber, and M. Woerner for 
detailed discussions and long-standing cooperation on the 
transport theory of QCLs. This work was supported by 
the Swedish Research Council (VR). 
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